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SECTION  1 


INTRODUCTION 


The  primary  objective  of  this  effort  has  been  the  implementation  of  a  non-reflecting  bound¬ 
ary  for  use  with  existing  finite-element  codes  to  perform  nonlinear  ground-shock  analyses  of 
buried  structures.  This  boundary  is  based  on  the  first-order  doubly  asymptotic  approxima¬ 
tion  (DAA\)  for  elastodynamic  scattering  f Geers  and  Yen  (1981),  Underwood  and  Geers 
(1981 )  .  In  addition,  a  staggered  solution  procedure  is  utilized  to  partition  the  global  equa¬ 
tions  in  order  to  achieve  both  computational  efficiency  and  software  modularity  [Felippa 
and  Park  (1980),. 

This  work  extends  that  of  Underwood  and  Geers  (1981)  for  linear  ground-shock  problems, 
wherein  the  DAA  surface  is  placed  on  the  surface  of  the  buried  structure.  Here,  the 
DAA  surface  is  moved  some  distance  out  from  the  surface  of  the  structure,  enclosing  both 
the  structure  and  a  portion  of  the  surrounding  soil  medium,  which  may  be  treated  with 
nonlinear  finite  elements.  Other  extensions  include  formulation  and  implementation  for 
general  2-D  and  3-D  problems,  improved  discretization  of  the  DAA  surface  with  higher- 
order  interpolation  functions,  and  utilization  of  a  conditionally  stable  staggered-solution 
procedure. 

1.1  Doubly  Asymptotic  vs.  Singly  Asymptotic  Approximations 

It  is  important  to  differentiate  between  doubly  asymptotic  approximations,  which  address 
quasi-static  and  wave-propagation  effects  simultaneously,  and  singly  asymptotic  approxi¬ 
mations,  which  address  these  effects  separately  [see,  e.g.,  various  papers  in  Kalinowski.  ed. 
(1981)  and  Datta.  ed..  (1982),  and  Cohen  and  Jennings  (1983)}  For  example,  representa¬ 
tion  of  the  external  medium  by  an  elastic  foundation,  which  may  be  quite  satisfactory  at 
low  frequencies,  does  not  account,  at  higher  frequencies,  for  energy  dissipation  through 
outward  propagation  of  scattered  waves.  On  the  other  hand,  representation  of  the  exter¬ 
nal  medium  by  a  viscous  boundary,  which  may  be  quite  satisfactory  for  wave-propagation 
problems,  does  not  provide  elastic  restoring  forces  in  the  static  limit. 

A  response-averaging  method  originally  proposed  by  Smith  (1974)  and  extended  by  Cun- 
dall.  et  a/.  (1978)  also  fails  in  the  static  limit.  For  example,  consider  the  response  of  a 
rigid  structure  surrounded  by  an  infinite,  linear-elastic  medium  to  an  internal,  quasi-static 
point  force.  A  computational  model  for  this  problem  might  consist  of  the  rigid  structure 
surrounded  by  a  portion  of  the  medium  enclosed  by  a  non-reflecting  boundary.  If  this 
boundary  is  that  of  Smith,  the  total  response  of  the  structure  is  the  average  of  two  re¬ 
sponses,  one  dependent  on  the  stiffness  of  the  bounded  portion  of  medium  enclosed  by 
a  rigid  boundary,  and  the  other  associated  with  the  structure  and  bounded  portion  of 
medium  floating  freely  in  space.  Unfortunately,  the  latter  response  grows  indefinitely  in 
the  static  limit  because  the  freely  floating  system  is  not  in  static  equilibrium.  In  contrast, 
doubly  asymptotic  approximations  approach  exactness  in  the  static  limit. 


1.2  Outline  of  Remainder  of  Report 

Section  2  of  this  report  derives  the  first-order  doubly  asymptotic  equations  of  motion  for 
a  buried  structure  excited  by  a  transient  incident  wave.  Section  3  deals  with  formation 
of  the  medium  stiffness  matrix  required  for  the  low-frequency  component  of  DAAi.  The 
staggered-solution  procedure  and  associated  stability  analysis  are  discussed  in  Section  4, 
which  establishes  the  time-increment  limitation  of  the  conditionally  stable  algorithm.  Sec¬ 
tion  5  describes  the  implementation  of  the  formulation  as  computer  software,  and  presents 
numerical  results  for  two  canonical  problems,  viz.,  excitation  of  an  infinite-cylindrical  aud 
a  spherical  shell  by  a  plane  dilatational  wave.  Section  6  concludes  the  report  with  some 
observations  and  recommendations  for  future  work. 


SECTION  2 


GOVERNING  EQUATIONS 


This  section  presents  the  governing  equations  for  the  finite-element  (FE)  model  of  the  struc¬ 
ture  along  with  a  portion  of  the  surrounding  soil  medium,  and  for  the  boundary-element 
model  (BE)  of  the  non-reflecting  DAA  surface.  These  equations  are  then  partitioned,  and 
a  staggered-solution  procedure  is  introduced  to  solve  for  transient  response.  Throughout 
the  development,  the  dependence  of  excitation  and  response  quantities  on  time  is  implicit. 


2.1  Finite-Element/Boundary-Element  Model 

Let  x  be  the  computational  vector  of  displacement  response  in  global  coordinates  for  the 
FE  model  of  the  structure  and  portion  of  surrounding  medium.  The  governing  equations 
for  the  finite-element  model  are  then  [see,  e.g.,  Zienkiewicz  (1977)} 


Msx  -f  Dsx  -  Ksx  =  (e  -i-  f, 


(2.1) 


where  Ms,  Ds  and  Ks  are  the  mass,  damping  and  stiffness  matrices,  respectively,  for 
the  FE  model.  fe  is  the  computational  vector  of  external  medium  forces  imposed  by  the 
DAA  surface  and  f,  is  the  vector  of  internal  nonlinear  forces:  as  usual,  a  dot  denotes 
differentiation  in  time.  Compatibility  of  forces  and  displacements  at  the  DAA  surface  may 
be  expressed  as  Geers  and  Underwood  (1981) 


fe  =  -Gg 
u  =  Grx 


(2.2) 


where  g  and  u  are  the  global  force  and  displacement  vectors,  respectively,  for  the  BE  model 
of  the  DAA  surface  and  G  is  the  force-transformation  matrix  from  BE  to  FE  coordinates. 

Now  the  force  vector  g  and  displacement  vector  u  may  be  decomposed  into  incident-wave 
and  scattered-wave  components  as 


S  g/  +  S5  (2.3) 

U  =  U/  —  Uy 

where  g,  is  the  known  force  vector  associated  with  a  free-field  incident  wave  and  g*  is 
the  unknown  force  vector  associated  with  the  wave  scattered  by  the  structure.  It  is  worth 
noting  that  this  dual  decomposition  does  not  require  constitutive  linearity  of  the  medium 
to  be  valid,  for  g<-  and  U5  may  each  be  viewed  as  merely  the  difference  between  two  vectors, 
one  obtaining  with  the  structure  absent  and  the  other  obtaining  with  the  structure  present. 

2.2  Doubly  Asymptotic  Approximation 

A  first-order  DAA  is  used  here  to  relate  the  scattered-force  vector  and  the  scattered- 
displacement  vector  U5  Geers  and  Yen  (1981)  and  Underwood  and  Geers  (1981)  .  This 


approximation  approaches  exactness  in  both  the  high-  and  low-frequency  limits,  and  effects 
a  smooth  transition  between.  The  development  of  DAAi  for  a  linear,  isotropic  external 
medium  proceeds  as  follows. 

At  high  frequencies ,  the  geometrical  vector  of  scattered-wave  surface  tractions  for  the  DAA 
surface  corresponding  to  normal  and  tangential  motions  of  that  surface  is  given  by 

t'sip)  =  PmCmu's{p)  (2-4) 

where  p  denotes  a  point  on  the  surface,  pm  is  the  mass  density  of  the  medium,  and  Cm  is 
the  diagonal  sound-speed  matrix  corresponding  to  u's ,  which  is  the  geometrical  vector  of 
normal  and  tangential  sc  altered- wave  velocities.  For  the  component  of  u's  normal  to  the 
DAA  surface,  the  corresponding  matrix  component  is  the  dilatational  velocity,  while  for 
each  component  of  u's  tangential  to  the  DAA  surface,  the  corresponding  matrix  component 
is  the  shear  velocity. 

Now  the  local-coordinate  vectors  of  (2.4)  may  be  transformed  into  global-coordinate  vectors 
as 

u'5(p)  =  Q(p)  u s(p),  t's(p)  =  Q(p)  ts{p )  (2.5) 

to  obtain,  inasmuch  as  Q_I  =  Ql .  where  the  superscripts  -1  and  t  denote  inverse  and 
transpose,  respectively. 

ts(p)  =  Ql(p)  PmCm  Q{p)  «s(p)  (2.6) 

Hence  boundary-element  discretization  of  us  as  see,  e.g..  Zienkiewicz  (1977) 

«s(p)  =  N(p)us  (2.7) 

where  N(p)  is  a  matrix  of  shape-functions  and  Us  is  a  vector  of  displacement  degrees  of 
freedom,  and  definition  of  the  high-frequency  scattered-wave  force  vector  as 

g s=  f  N‘(p)  ts(p)  dS  (2.8) 

yield,  for  high-frequency  motions. 

g  s  =  Dmu.s  (2.9) 

in  which 

Dm  =  f  N1  Q<  Prr,CmQ-NdS  (2.10) 

At  low  frequencies,  the  scattered-wave  force  computational  vector  is  given  by  the  quasi¬ 
static  relation 


g  1S  =  Kmus 


(2.11) 


where  Km  is  a  full,  nonsymmetric  stiffness  matrix  for  the  boundary-element  mesh,  whose 
construction  is  described  in  the  next  section. 


Finally,  the  first-order  doubly  asymptotic  approximation  DAA/  is  formed  by  the  superpo¬ 
sition  of  gls  and  g$  to  obtain 


gs  =  Dmu5  -  Kmus 


(2.12) 


Now  the  assumption  embodied  in  DAA!  of  a  constitutively  linear  medium  for  the  scattered 
wave  is  justified  within  the  framework  of  classical  plasticity  theory  if  the  material  point 
for  every  exterior  location,  «.e.,  every  location  in  the  medium  outside  the  DAA  surface, 
remains  within  its  corresponding  yield  surface  when  and  after  the  scattered  wave  arrives 
at  the  DAA  surface.  For  incident  waves  with  sufficiently  rapid  decay  rates  and  for  a  DAA 
surface  sufficiently  removed  from  the  surface  of  the  structure,  the  scattered  wave  causes 
minor  perturbations  about  an  elastic  state  at  each  exterior  location,  thereby  satisfying  the 
preceding  condition. 


The  assumption  of  material  isotropy  outside  the  DAA  surface  cannot  be  rigorously  main¬ 
tained  if  the  material  has  suffered  plastic  excursions  in  response  to  the  incident  wave. 
However,  it  is  likely  that  the  resulting  anisotropy  is  no  more  pronounced  than  that  char¬ 
acterizing  the  ambient  state,  which  is  generally  uncertain  in  practical  cases.  Hence,  while 
an  extension  to  material  orthotropy  may  be  theoretically  possible,  it  may  not  be  worth 
the  trouble. 


2.3  Response  Equations 


Introduction  of  the  first  of  (2.2)  and  (2.3)  into  (2.1)  and  of  the  second  of  (2.2)  and  (2.3) 
into  (2.12)  yields  the  doubly  asymptotic  equations  of  motion 


M?x  -  Dfx  -  K,x  =  -G{g,  -  gs}  -  f, 


(2.13) 


gs  =  Dm{Grx  -  in}  -  Km{Grx  -  u/} 


which  may  be  numerically  integrated  in  time  to  obtain  the  solution  vectors  x  and  gs. 
Because  M?,  and  K..,  are  typically  large  and  banded,  while  Km  is  relatively  small  and 
full,  it  is  not  computationally  practical  to  introduce  the  second  of  these  equations  into  the 
first  to  eliminate  g5. 


However,  because  Dm  is  banded  and  multiplies  the  highest -derivative  terms  in  the  second 
of  (2.13).  it  is  advantageous  to  apply  the  technique  of  augmentation  Park,  or  a/.  (1977)  . 
which  here  merely  involves  introducing  the  second  of  12.13)  into  the  hrst.  moving  the  term 
containing  to  the  left  side  of  the  resulting  set  of  equations,  and  keeping  G  K„,G  x  on 
the  right.  This  yields  the  augmented  doubly  asymptotic  equations  of  motion 
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SECTION  3 


!  MEDIUM  STIFFNESS  MATRIX 

» 

» 

\  This  section  describes  the  construction  of  the  boundary-element  stiffness  matrix  that  re- 

I  lates  the  scattered-wave  force  and  displacement  vectors  at  low  frequencies.  The  develop- 

I  ment  is  based  on  Somigliana’s  identities,  which  derive  from  Betti’s  reciprocal  work  theo¬ 

rems  and  Kelvin's  problem  of  a  point  load  in  an  infinite  elastic  medium  [see,  e.g.,  Kupradze 
(1965),  Rizzo  (1967),  Cruse  (1969),  Lachat  and  Watson  (1976)] 

3.1  Elastostatic  Boundary-Integral  Equations 

The  surface  behavior  of  an  elastic  medium,  whether  occupying  an  exterior  or  interior 
region,  may  be  expressed  as  [Rizzo  (1967),  Cruse  (1969)] 


c(p)u(p) -T  J  T(p,  q)u(q)  dT  q  =  J  U(p,  q)t  (q)  dT  q  (3.1) 

where  p  is  a  point  on  the  boundary  and  q  is  the  integration  variable,  and  where  u(p)  and 
t(p)  are  d  x  1  vectors  (d  =  2  or  3)  of  medium  displacements  and  tractions  in  Cartesian 
coordinates  on  the  boundary  at  p.  The  elements  Tl}(p,  q)  and  f7t;(p,  q)  of  the  d  x  d  matrices 
T(p,q )  and  U(p,q)  are  fundamental  solutions  for  the  tractions  and  displacements  at  a 
location  q  in  the  direction  i  due  to  a  point  load  at  location  p  in  direction  j.  With  <5t;  as 
the  Kronecker  symbol,  each  element  of  the  matrix  c  is  defined  as 


c.j(p)  =  kSi}  (3.2) 

if  there  exists  a  continuous  tangent  at  p.  or.  with  T(  as  the  surface  of  a  sphere  of  radius  e 
centered  at  p, 


c.j(p) 


(3.3) 


if  the  tangent  is  not  continuous.  A  simple  method  for  the  evaluation  of  c,,  is  given  in 
Appendix  A. 

Now  an  element  of  the  two-dimensional  displacement-kernel  matrix  i'(p.q)  for  plane-strain 
problems  is  given  by 


Ul}(p,q) 


-1 

S7r(]~-~u)C 


,(3 


4  i/)/n(r)6,; 


(3.4) 


where  <J  and  v  are  the  shear  modulus  and  Poisson’s  ratio. respectively,  and  r  —  r(p.q)  is 
the  distance  between  the  load  point  p  and  the  field  point  q:  the  derivatives  are  taken  with 
reference  to  the  coordinates  of  q.  With  p,  and  q,  as  t  he  coordinates  of  p  and  q.  respectively. 


D  =  9,  -  P, 


t  =  (r,rt)  2 
9»  -  Pt 


In  contrast,  an  element  of  the  three-dimensional  displacement-kernel  matrix  U(p,q)  is  given 
by 

Utj{p,q)  =  ;<3  “  At/)6'J  ^  r ’«  rv  1  (3-6) 

Finally,  an  element  of  the  traction-kernel  matrix  T(p,  q)  for  both  two-  and  three-dimensional 
problems  is  given  by 


Ttj{p,q)  =  4ct7^  _  ^^{K1  “  2u)6'}  +  Pr"  r'i  )r’*n*  ~  (!  ~  2i/)(r„ny  -  r,;  nt)}  (3.7) 

where  a,  and  n;  are  direction-cosines  for  the  surface  normal  at  q.  The  two-  and  three- 
dimensional  forms  are  explicitly  obtained  by  letting  a  =  1,  2  and  0  =  2,  3,  respectively. 

3.2  Discretization 

Numerical  solution  of  the  integral  equation  (3.1)  requires  discretization  of  the  DAA  surface, 
over  each  boundary  element  of  which  the  displacement  and  traction  vectors  are  approxi¬ 
mated.  The  curved  isoparametric  elements  of  finite-element  theory  offer  both  the  generality 
and  the  accuracy  needed  for  this  purpose.  With  this  approach,  the  global  Cartesian  coor¬ 
dinates  of  any  point  in  an  element  are  taken  as  related  to  the  nodal  coordinates  by  c.f. 
(2.7) 

x[p)  =  N(p)  x  (3.8) 

i.c.,  the  same  shape  functions  are  used  to  approximate  element  geometry,  displacements 
and  tractions.  This  allows  interpolated  displacements  and  tractions  along  the  DAA  curve 
in  two-dimensional  space  to  be  integrated  over  a  normalized  length  in  ^-coordinate  space, 
and  similar  quantities  over  the  DAA  surface'  in  t hroe-dimensional  space  to  be  integrated 
over  a  standard  2  *  2  normalized  square  in  £i.  ^-.--coordinate  space. 

On  an  element-by-element  basis,  (3.8)  becomes 


xe(f)  =  xe* 


where  re(£'’)  is  the  d  ■  1  vector  of  Cartesian  coordinates  of  a  point  in  element  e.  the  A 'k(C) 
are  the  element  shape  functions,  and  x'k  is  the  d  ■  1  vector  of  Cartesian  coordinates  of  the 
A- 1 h  element  node:  also.  ^  in  2-D.  but  £  in  3-D.  The  elements  used  in  this 


study  are  the  three-noded,  quadratic,  curved  element  for  2-D  analysis  and  the  eight-noded, 
quadratic,  serendipity  element  for  3-D  analysis.  The  shape  functions  for  the  three-noded 
quadratic  element  are 


I 

ij 


=  1  -  f2 

iv3=  |e(€  +  i) 


(3.10) 


where  £  €  [-1,1);  the  nodes  are  located  at  f  =  -1,0,1.  The  shape  functions  for  the  eight- 
noded  quadratic  element  are 

Ni  =  -i(i  -  ei)(i  -  €2)(i  +  6  +  &) 

^2  -  |(1  -  tf)(l  -  &) 

^3  =  ^(1  +  £l)(l  -  f2)Ul  -  ^2  ~  1) 

n4  =  5(1  +  €i)(i  -  €2)  /„  in 

^5  =  4(1  +  Ci)(i  -  £2)(£i  +  £2  -  1) 

.v6  =  |(i  -  )(i  -  6) 

i(l-(,)(lt(2)K,-s£2-  1) 

7V8  =  I(l-6)(l-d) 


where  €  1-1,1  and  £2  €  j-1,1  j,  and  all  nodes  lie  at  the  intersections  of  the 
and  the  £2  =  -1,0,1  lines,  except  at  0,0.  where  there  is  no  node. 

3.3  Matrix  Assembly 

With  DAA-surface  coordinates,  displacements  and  tractions  approximated  as 

x(p)=N(p)x  u(p)  =  N(p)u  t(p)  =  N(p)t 
(3.1)  may  be  expressed  at  a  node  P  as 

E  f 

c{P)u(P)-Y  T(P.qC')  V.V,(n<  J(C)dC 
, _ ,  >r 


=  -1.0.1 


(3.12) 


-  V 


f  l'[P.q  r  )  y.\k{?)t'kJ(£)dr 


(3.13) 


where  E  is  the  total  number  of  elements  on  the  DAA  surface  and  ,/(£e)  is  the  Jacobian 
for  xe  -  transformation:  also.  d(f  -  d£e  in  2-D.  but  d£  =  d^d^n  in  3-D.  Finally, 
coalescence  of  element  contributions  at  common  nodes  is  implicit  in  (3. 1 3).  The  numerical 
techniques  used  to  evaluate  the  integral'-  in  this  equation  are  discussed  in  Appendix  A. 

Evaluation  of  (3.13)  at  every  node  on  the  DAA  surface  \  ields  a  se>  of  simultaneous  algebraic 
equations  that  can  be  expressed  in  the  form 
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Now  the  nodal  force  vector  g  corresponding  to  a  traction  distribution  t  on  the  DAA  surface 
is  given  by 

g  =  !  N‘(p)  t(p)  dV  (3.16) 

Introduction  of  the  third  of  (3.12)  and  of  (3.15)  into  this  relation  then  yields 

g  =  Kmu  (3.17) 

where  the  generally  non-symmetric  medium  stiffness  matrix  Km  is  given  by 

Km  =  [  J  N4N  dr  B  *A  (3.18) 

A  symmetric  form  may  be  obtained  as 

Km  =  i(Km-K‘j  (3.19) 

which  is  identical  to  that  derived  from  energy  considerations  ' Zienkiewicz.  Kelly  and  Beti- 
ess  (1977) j.  As  indicated  in  Appendix  B.  however,  the  use  of  Km  generally  yields  numerical 
results  inferior  to  those  produced  bv  Km- 


SECTION  4 


STAGGERED  SOLUTION  PROCEDURE 


In  the  interest  of  computational  efficiency,  the  augmented  doubly  asymptotic  equations 
of  motion  given  by  (2.14)  are  solved  with  a  staggered  solution  procedure.  The  procedure 
is  conditionally  stable,  requiring  that  the  time  increment  be  smaller  than  the  shortest 
medium- boundary  period  divided  by  n.  This  shortest  period  may  be  obtained  by  determin¬ 
ing  the  highest  natural  frequency  for  the  eigenproblem 

w2Mfx=GKmG!x  (4.1) 

In  cases  where  the  surrounding  soil  does  not  appreciably  stiffen  the  embedded  structure  be¬ 
yond  its  inherent  level,  the  highest  medium-boundary  frequency  is  substantially  lower  than 
the  highest  natural  frequency  characterizing  the  structure  itself,  thereby  allowing  the  ana¬ 
lyst  to  carry  out  stable  calculations  with  a  relatively  large  time  increment.  The  remainder 
of  this  section  describes  the  staggered-solution  procedure  and  the  stability  analysis  that 
leads  to  (4.1). 

4.1  Solution  Algorithm 

To  construct  the  staggered  solution  procedure  for  (2.14).  those  equations  are  expressed  at 
mid-step  as 


2  -  DTxn  +  1/2  -  K.oXn-,  2  =  f„_  l '2  ~  KM  xn_1/2  (4.2) 

where  the  time  step  n  =  t/  At.  in  which  t  and  At  are  time  and  fixed  time  increment,  re¬ 
spectively,  and  where  the  total  damping  matrix  Dj-.  the  medium-boundary  stiffness  matrix 
Km,  and  the  total  force  vector  f  are  given  by 

Dr  =  Ds  ■*-  GDmG' 

Km  =  GKmG‘  (4.3) 

f  -  -Gg/  —  GDr^ii/  —  G  Kmu/  f, 

The  integration  algorithm  utilized  is  the  trapezoidal  rule  see.  e.g..  Henrici  (1962)  .  for 
which 


Xn-t  2  ~  (xr,  -  1 

2  -  xn)  6 

*«-l  2  =  (X„- i 

2  -  X„)  b 

Xn- 1  =  2xn_ ] 

2  ~  Xn 

x n —  1  2xn  _ ] 

2  Xn 

(4.4) 


where  b  —  At  2.  Introduction  of  the  first  and  then  the  last  of  these  into  the  third  yields 
the  standard  form 


1 1 


Now  the  first  two  of  (4.4)  are  introduced  into  the  left  side  of  (4.2)  and  xn+lj/2  on  the  right 
side  of  (4.2)  is  predicted  as  xrn^x,n  to  obtain  the  set  of  algebraic  equations 


E,xn.1/2  -  en+1/2  -Erx^1/2 

where 

Ei  =  Ms  -  <5Dr  -  <52K^ 

Er  =  62  Km 

Gn+l/2  =  <52fn^l/2  ~  Ms(x„  +  <5xn)  -  <5 DrXn 
Finally,  the  prediction  x^_1/0  is  based  on  the  one-term  extrapolation 

Xn-l/2  =  X” 


(4.6) 


(4.7) 


(4.8) 


The  preceding  staggered  solution  procedure  leads  to  the  following  computational  sequence 
to  determine  system  response  at  time  step  n  —  1: 


(a) 

f"n-M/2  “  (fn  ~  fn-l)/2 

(b) 

e„.i/2  =  ^ 2 -f  i /' 2  -  Ms(xn  -  6xn) 

-  6  DTX„ 

(c) 

—  X 

1/2  n 

(d) 

xn~l/2  —  E(  1/2  —  ErXni|  ,, 

(e) 

Xn-  1  —  2xn4.i  '2  ^ n 

(!) 

Xn  —  1/2  =  (xn  — 1/2  -  xn)  '  b 

(9) 

xn- 1/2  =  1(fn-^l/2  _ 

-  Krx„^i/2 

(h) 

xn-l/2  =  xn  ~  <5xn*l  2 

(i) 

X„_  i  -  2x„^  i  2  -  X„ 

where  the  total  stiffness  matrix  Kf  =  K,.  -  K.vj-  To  improve  accuracy,  an  iterative  loop 
has  been  introduced  at  (dj.  wherein  x;n_12  on  the  right  is  corrected  to  the  previously 
calculated  value  of  x„  +  1/2;  two  iterations  generally  produce  satisfactory  convergence.  The 
calculation  starts  at  n  —  0  with  x,,  =  x,,  =  0. 

4.2  Stability  Analysis 

Park  (1980)  has  performed  a  stability  analysis  of  a  generalized  form  of  the  staggered 
solution  procedure  just  described.  The  result  is  that  the  procedure  is  computationally 
stable  if  no  root  of  the  characteristic  equation 


v  \  "  v ;  '■ 

•/i'  -  /.  •. 


det  ]  z2(Mf  -  62Km)  ~  zSDt  ~  d‘Kr  0  (1.9) 

has  a  positive  real  part.  Verification  of  this  condition  is  relatively  straightforward  when  all 
of  the  matrices  in  (4.9)  are  symmetric;  it  is  generally  quite  difficult  when  one  or  more  is  not . 
Unfortunately,  as  discussed  in  Section  3,  the  medium  stiffness  matrix  Km  is  nonsymmetric. 
which  pollutes  Km  and  K7.  Fortunately,  however.  Km  constitutes  a  small  perturbation 
of  Km,  which  is  symmetric:  hence  it  is  appropriate  to  consider  the  characteristic  equation 

det  1  z  ~  ( M  ^  =0  (4.10) 

where  Km  =  G  KmG(  and  K7-  =  Kf  -  Km- 

As  discussed  on  page  255  of  Bellman  (1970).  no  root  of  (4.10)  has  a  positive  real  part  if 
(M*  -  62Km)<  Dt  and  K7-  are  all  non-negative  definite  and  either  (Ms  -  62Km)  orKr 
is  positive  definite.  On  physical  grounds.  Dr  and  K7  are  both  non-negative  definite,  but 
generally  not  positive  definite.  However,  inasmuch  as  Mf  is  positive  definite,  (M.5  — <52Km) 
is  positive  definite  if  6  is  sufficiently  small.  The  degree  of  smallness  defines  the  stability 
requirement,  as  discussed  next. 


Consider  the  following  first  etgenproblem: 


Qx  =  A  x 


(4.11) 


where  Q  =  Ms  'k^.  This  problem  yields  non-negative  real  eigenvalues  and  real  eigen¬ 
vectors.  These  eigenvectors  may  be  assembled  into  a  modal  transformation  matrix  ¥  that 

_  *  .  *  *  d  1 

diagonalizes  Q  as  *  Q*  =  Q  and  normalizes  as  9  =  I,  the  identity  matrix.  Hence 
the  introduction  into  (4.11)  of  a  transformation  from  physical  to  generalized  coordinates 
as  x  =  ^  y  and  subsequent  premultiplication  through  by  ft*  yield  the  diagonal  eigenvalue 


matrix 


Ao  -  Q 


Consider  next  the  following  second  eigenproblem: 


K  AM  x 


(4.12) 


(4.13) 


whose  eigenvalues  and  eigenvectors  are  the  same  as  those  of  the  first  eigenproblem.  Hence 
the  transformation  from  physical  to  generalized  coordinates  and  premultiplication  through 
bv  vields 


A  km  =  (M‘)-’Km 

•  -  d  ,  ■ 

where  M  M  AP  and  K  X1  *  ¥  K  vr®  A  k  V/  is-  of  course,  identical  1o  A  q. 

Finally,  consider  t he  following  third  eigenproblem : 


(A14) 


' .'.'A,' >  ......  •  .  . 

‘/‘cVv'  /  . 

>  V'  -  ' 

•  *  v  -V.'  •  n.  «.  -  "1  • 


■  *  -*  ;w.« 


/  s  f 


(Ms  -  62Km)  x  =  Ax 

Transformation  and  premultication  through  as  before  yields 

AM_K  =  Mds-62KM 

=  Mfl-62AK/M 

=  Ml!  -  62An 


(4.15) 


(4.16) 


Hence  the  eigenvalues  of  (M,  -  62 K^)  are  all  positive,  and  thus  (Mf  -  <52K«)  is  positive 
definite,  if  6 2  times  the  largest  eigenvalue  A qax  is  less  than  unity.  With  A QaT  —  (t jQax)2. 
this  yields  the  stabthty  requirement 


which  is  stated  in  slightly  different  terms  at  the  beginning  of  this  section. 

Establishment  of  the  stability  requirement  (4.17)  for  a  symmetric  medium  stiffness  matrix 
facilitates  the  estimation  of  a  similar  requirement  for  a  non-svmmetric  one.  Clearly,  no  root 
of  (4  9)  has  a  positive  real  part  if  6  is  vanishingly  small,  as  M.,  is  symmetric  and  positive 
definite,  and  is  symmetric  and  non-negative  definite.  Also,  on  physical  grounds,  the 
eigenvalues  of  (Ms)-1Km  must  be  real  and  non-negative.  Finally,  the  eigenvalues  for  the 
three  eigenproblems  above  differ  only  slightly  from  their  counterparts  when  K m  is  replaced 
by  K«  because,  as  illustrated  in  Appendix  B,  Kw  constitutes  a  small  perturbation  ofKm. 
Hence,  as  6  is  increased  from  zero,  all  the  roots  of  (4.9)  contain  negative  real  parts  until 
the  stability  requirement  (4.17)  is  approached,  where  u ,'qoz  now  pertains  to  the  use  ofKm. 
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SECTION  5 

IMPLEMENTATION  AND  COMPUTATION 

This  section  describes  the  techniques  used  to  implement  in  software  the  approach  delin¬ 
eated  above,  and  presents  numerical  results  generated  by  that  software.  Modern  software- 
engineering  techniques  are  used  jfeiippa  (1981)].  in  order  to  facilitate  extension  to  large- 
scale  production  analysis.  The  numerical  results  pertain  to  canonical  problems  involving 
plane,  dilatational  step-waves  that  envelop  infinite-cylindrical  and  spherical  shells  (Figure 
1).  These  problems  possess  known  analytical  solutions. 

5.1  Software  Implementation 

The  approach  described  in  Sections  2,  3  and  4  is  embodied  in  an  assembly  of  of  four 
software  entities: 

1.  Structural  Matrix  Generator:  The  structural  mass  and  stiffness  matrices.  M?  and  Kf 
in  (2.14),  are  generated  by  the  finite-element  code  DIAL  Ferguson  and  Cyr  (1984),: 
Ds  is  neglected.  The  structural  matrices  and  related  data  are  read  into  a  NICE  global 
database  ; Felippa  (1982)). 

2.  Medium  Matrix  Generator:  The  medium  damping  and  stiffness  matrices.  Dm  and 
Km  in  (2.14),  are  generated  by  software  developed  as  part  of  tHs  study  in  the  manner 
described  above:  the  force-transformation  matrix  G  is  constructed  as  a  correspondence 
table.  These  data  are  read  into  the  NICE  global  database. 

3.  Incident  Field  Generator:  The  incident-wave  displacement,  velocity  and  force  vectors. 
U/.  U/  and  g;  in  (2.14),  are  also  generated  by  software  developed  as  part  of  this 
study  in  the  manner  described  below:  as  these  are  time-dependent  vectors,  they  are 
calculated  dynamically  as  the  calculation  proceeds,  f,  is  taken  as  zero. 

4.  Staggered  Solution  Procedure:  The  solution  algorithm  described  in  Subsection  4.1  is 
implemented  as  a  NICE  procedure  using  a  command  language  interpreter  Felippa 
(1983)  .  The  matrix  operations  embedded  in  the  algorithm  are  performed  with  a 
matrix  utility  processor  for  data  in  unblocked  skyline  format  [Felippa  (1978).. 

The  FE  and  BE  models  are  constructed  independently,  although  the  element  grids  match 
at  the  boundary.  Geometrical  symmetry  is  exploited  in  both  canonical  problems. 

5.2  Incident-Wave  Vectors 

A  plane,  dilational  step-wave  characterized  by  a  velocity  jump  V’t,  and  propagating  in  the 
indirection  may  be  described  in  terms  of  a  scalar  potential  as 


<t>‘  =  ~ 


x i  -  a)'  H(cdt  -  X]  -  a) 


where  c,j  is  the  dilatational  speed  in  the  elastic  medium.  H  is  the  Heaviside  operator,  and 
-a  is  the  point  on  the  inaxis  where  the  wave  front  is  located  at  t  -  0.  The  applica¬ 
tion  of  classical  continuum  formulas  Achenbach  (1973)  yields  for  the  components  of  the 
geometrical  displacement  and  velocity  vectors  for  the  incident  wave 


H 


it** 

sss 


.N.-.  /■  u  .N.N.N.N  A.vGVaG'A''  ■/aV-.V  tfV.N.1 


e  F.Vfr>j 


(5.2) 


uj  =  <5lt  —  (cat  -  xx  -  a)H(cdt  -  xi  -  a) 

Cd 

u\  =  6uv0  H{cdt  -  X]  -  a) 

Hence  the  elements  of  the  computational  vectors  U/  and  U/  are  given  by  (6.2)  evaluated 
at  the  surface  nodes. 

Similarly,  the  components  of  the  incident-wave  stress  tensor  and  geometrical  surface- 
traction  vector  are  given  by  [Achenbach  (1973) j 

y 

o{.  =  -6tJ—  (A  +  2p6u)  H[cdt  ~  X]  -  a)  (5.3) 

Cd 

tj  = 

where  A  and  p  are  the  Lame  constants  and  the  n,  are  the  direction-cosines  for  the  surface 
normal.  Hence  the  computational  vector  g,  is  given  by  (3.16). 

5.3  Infinite  Cylindrical  Shell 

The  first  canonical  problem  is  that  of  an  infinite  cylindrical  shell  embedded  in  an  elastic 
medium  and  excited  by  a  transverse,  plane,  dilatational  wave  Garnet  and  Crouzet-Pascal 
(1966)].  The  parameter  ratios  for  this  problem  are  E^'Em  =  2.5  (Young's  modulus). 
h/a  =  0.01  (shell  thickness-to-radius),  PYiPm  -  1156  (mass  density).  vm  -  0.25  and 
I/.,  =  0.2  (Poisson’s):  these  pertain  to  a  concrete  shell  in  slow  granite.  The  duration  of  the 
rectangular  incident-wave  pulse  is  Cdt/a  =  10.  A  curved,  three-noded  shell  element  is  used 
to  model  the  shell,  so  that  the  FE/BE  discretization  employs  conforming  elements. 

The  first  computational  model  for  this  problem  places  the  DAA  boundary  directly  on  the 
shell  in  the  manner  of  Underwood  and  Geers  (1981).  The  use  of  six  curved  quadratic  ele¬ 
ments  over  the  half-model  yields  results  that  are  virtually  identical  to  those  of  Underwood 
and  Geers  (1981).  which  were  generated  with  twenty  linear  elements  over  the  half-model. 
Figure  2  shows  DAA  and  exact  displacement-response  histories;  agreement  is  seen  to  be 
excellent. 

The  second  computational  model  introduces  eight-noded  medium  finite  elements  between 
the  shell  and  the  DAA  boundary,  which  is  located  one  shell  radius  out  from  the  shell  surface 
(Figure  3).  The  displacement-response  histories  thus  produced  are  shown  in  Figure  -1  as 
solid  lines,  along  with  their  DAA  counterparts  from  Figure  3.  which  are  shown  as  dashed 
lines.  It  is  seen  that  the  use  of  medium  finite  elements  degrades  solution  accurac>  somew'hat 
by  introducing  spurious  oscillations  caused  by  ringing  of  the  mesh.  A  third  computational 
model,  which  locates  the  DAA  boundary  three  shell  radii  out  from  the  shell  surface,  yields 
results  that  are  even  more  oscillatory,  although  peak-response  values  are  still  satisfactory. 

5.4  Spherical  Shell 

The  second  canonical  problem  is  that  of  a  spherical  shell  embedded  in  an  elastic  medium 
and  excited  by  a  plane  dilatational  wave  Grafton  and  Fox  (1965).  Geers  and  Yen  (1981)  . 
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The  parameter  ratios  for  this  problem  are  the  same  as  those  for  the  infinite  cylindrical 
shell,  and  the  duration  of  the  rectangular  incident-wave  pulse  is  also  c<it  fa  =  10.  An  eight- 
noded  Ahmad  shell  element  is  used  to  model  the  shell,  so  that  this  FE/BE  discretization 
also  employs  conforming  elements. 

As  previously,  the  first  computational  model  for  this  problem  places  the  DAA  boundary 
directly  on  the  shell;  six  eight-noded  quadratic  elements  are  used  over  the  quarter-model  of 
the  shell  (Figure  5).  DAA-based  displacement-response  histories  are  compared  with  their 
exact  counterparts  in  Figure  6,  the  latter  having  been  generated  in  the  manner  of  Geers 
and  Yen  (1981).  Here  too,  agreement  is  seen  to  be  excellent. 

The  second  computational  model  introduces  twenty-noded  medium  finite  elements  between 
the  shell  and  the  DAA  boundary,  which  is  located  one  shell  radius  out  from  the  shell  surface 
(Figure  7).  The  displacement-response  histories  thus  produced  are  shown  in  Figure  8  as 
solid  lines,  along  with  their  DAA  counterparts  from  Figure  6.  which  are  shown  as  dashed 
lines.  Here  too,  it  is  seen  that  the  use  of  medium  finite  elements  degrades  solution  accuracy 
by  introducing  spurious  oscillations  caused  by  ringing  of  the  mesh. 
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SECTION  6 


CONCLUSION 

This  report  has  documented  the  formulation  and  implementation  of  a  non-reflecting  bound¬ 
ary  for  use  with  existing  finite-element  codes  to  perform  nonlinear  ground-shock  analyses 
of  buried  structures.  The  boundary  is  based  on  a  first-order  doubly  asymptotic  approxi¬ 
mation  (DAA\)  for  disturbances  propagating  outward  from  a  selected  portion  of  the  soil 
medium  surrounding  the  structure  of  interest.  The  resulting  set  of  first-order  ordinary 
differential  equations  is  then  combined  with  the  second-order  equations  of  motion  for  the 
finite-element  model  so  as  to  facilitate  solution  by  a  staggered  solution  procedure.  This 
procedure  is  shown  to  be  computationally  stable  as  long  as  the  time  increment  is  smaller 
than  a  limiting  value  based  on  the  finite-element  mass  matrix  and  the  DAA-boundarv  stiff¬ 
ness  matrix.  Computational  results  produced  by  the  boundary  are  compared  with  exact 
results  for  linear  canonical  problems  pertaining  to  infinite-cylindrical  and  spherical  shells. 

6.1  Observations 

It  is  appropriate  here  to  offer  some  comments  regarding  the  work  described  above: 

1.  As  pointed  out  in  the  Introduction,  doubly  asymptotic  approximations  are  clearly 
superior  to  singly  asymptotic  approximations,  the  former  incorporating  both  radiative 
energy  dissipation  and  elastic  restoring  forces,  the  latter  accounting  for  only  one  or 
the  other. 

2.  While  the  medium  damping  matrix  may  be  interpreted  in  terms  of  local  dashpots 
positioned  on  the  DAA  surface,  the  medium  stiffness  matrix  is  not  so  easily  regarded: 
attempts  to  simplify  the  fully  coupled  nature  of  Km  merely  degrade  the  validity  of 
the  low-frequency  approximation. 

3.  Although  it  is  tempting  to  use  a  symmetric  medium  stiffness  matrix  in  DAA  compu¬ 
tations,  the  resulting  loss  of  accuracy  constitutes  too  high  a  price. 

4.  The  computational  stability  requirement  (4.17)  is  a  generous  one  when  the  soil  is 
substantially  softer  than  the  structural  material:  when  this  is  not  the  case,  however, 
more  efficient  computations  might  be  realized  with  an  unconditionally  stable  staggered 
solution  procedure,  which  is  yet  to  be  developed. 

5.  The  use  of  modern  software-engineering  techniques,  as  embodied  in  the  SICE  In¬ 
tegrated  Software  System,  greatly  facilitates  the  implementation  of  methods  for  the 
analysis  of  coupled  systems. 

6.  The  results  for  the  linear  canonical  problems  once  again  demonstrate  the  difficulty  of 
propagating  a  discontinuous  wave  front  through  a  finite-element  grid  and.  in  contrast  , 
the  good  performance  of  a  boundary-element  grid  located  directly  on  the  surface  of 
the  structure. 

6.2  Future  Work 

Future  R&D  work  in  this  area  could  profitably  pursue  the  following  paths: 


The  usefulness  of  the  non-reflecting  DAA\  boundary  in  nonlinear  problems  should  be 
more  stringently  assessed  by  applying  it  to  nonlinear  canonical  problems;  the  challenge 
here  is  to  find  “exact”  solutions  for  such  problems  against  which  to  compare  the 
approximate  solutions. 

A  non-reflecting  DAA  boundary  should  be  developed  for  a  medium  half-space,  this 
for  application  in  near-surface  ground-shock  analyses. 

An  unconditionally  stable  staggered  solution  procedure  should  be  formulated  for  prob¬ 
lems  not  amenable  to  the  conditionally  stable  procedure. 

A  new  approach  should  be  sought  for  satisfactorily  propagating  discontinuous  wave 
fronts  through  finite-element  grids;  failing  this,  the  option,  in  nonlinear  response  prob¬ 
lems.  of  placing  the  DAA  grid  directly  on  the  surface  of  the  structure  !  Underwood 
and  Geers  (1980)]  should  be  revisited. 
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Figure  4.  Displacement  response  histories  for  the  infinite  cylinderical  shell 
(solid  curves:  finite  elements  around  shell: 
dashed  curves:  DAA  boundary  on  shell  surface). 
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Figure  7.  Quarter-model  grid  for  spherical  shell  problem 
(Finite  elements  extending  to  r  =  2a). 
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ERROR  IN  CONSISTENT  FORCES  GENERATED  ON  A  SPHERICAL  CAVITY 
BY  A  UNIFORM  IMPOSED  DISPLACEMENT 
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o  s  Non  symmetric 


Figure  9.  Error  in  nodal-force  values  produced  by  symmetric  and  nonsymmetric 

medium  stiffness  matrices. 


APPENDIX  A 


NUMERICAL  INTEGRATION  TECHNIQUES 

Discretization  of  the  DAA  boundary  makes  it  possible  to  approximate  (3.1)  by  a  system  of 
linear  algebraic  equations  for  nodal  values  of  surface  displacement  and  traction,  t.e.,  (3.14). 
The  coefficients  in  these  equations  are  obtained  by  integrating,  by  means  of  quadrature 
formulas,  products  of  kernel  functions  and  shape  functions  over  the  boundary  elements,  as 
indicated  in  (3.13).  In  this  regard,  it  is  necessary  to  distinguish  between  to  fundamentally 
different  types  of  integrals  that  arise. 

The  first  type  of  integral  occurs  when  the  node  P  does  not  belong  to  the  element  over 
which  the  integral  is  being  performed.  This  type  is  regular ,  because  the  integrand  varies 
smoothly  over  the  surface.  Simple  Gaussian  quadrature  formulas  may  then  be  used.  In 
two  dimensions. 


f  +  i  M 

/  /U)dsc~  £>/  ftti 

J  -1  /=i 


where  the  w\  are  weighting  factors,  the  £/  are  the  coordinates  of  the  integration  points  and 
M  is  the  total  number  of  integration  points.  Similarly,  in  three  dimensions. 
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The  second  type  of  integral  occurs  when  the  node  P  belongs  to  the  element  over  which  the 
integral  is  being  performed.  This  type  is  singular,  because  the  integrand  grows  without 
bound  at  P.  The  techniques  used  to  evaluate  the  singular  integrals  encountered  in  this 
study  are  described  below. 

A.l  Singular  Integrals  Involving  the  Traction  Kernel 

For  this  singular  case,  there  exists  no  quadrature  formula  suitable  for  the  calculation  of 
the  integral  T,r  The  coefficient  of  this  integral  for  the  singular  node  together  with  the  c!; 
term  form  the  leading  diagonal  submatrix  of  coefficients  of  u}  in  equation  (3.13).  These 
coefficients  can  be  expediently  calculated  by  noting  that  a  stress  field  corresponding  to  a 
rigid  body  translation  of  the  body  is  zero.  In  this  case  equation  (3.14)  becomes 

A  u  -  0  (.4.3) 

where  u  is  a  vector  of  unit  rigid  body  displacements.  The  diagonal  terms  of  A  are  simply- 
given  by 

o„  1  ^  a,,  (-4.4) 


A. 2  Singular  Integrals  Involving  the  Displacement  Kernel 

In  two  dimensions,  a  quadrature  rule  based  on  the  theory  of  [Takahasi  and  Mori  (1973)] 
is  utilized  integrate  the  ln(r)  singularity  contained  in  the  UtJ  kernel.  Such  a  quadrature 
rule  has  been  successfully  used  for  two-dimensional  acoustic  scattering  problems  |J3urton 
1976],  First,  let  the  integration  of  this  singularity  in  the  intrinsic  coordinate  system  be 
represented  by 

1  =  1  (A5) 

where  /(£)  may  have  singularities  at  ±1.  Then  the  value  of  this  integral  is  given  by  the 
following  quadrature  formula 

N 

7  a:  V  wnf(‘n)  (A. 6) 

n  =  -  A’ 
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The  values  of  N  =  4  and  ft  =  0.75  were  used  to  construct  a  9-point,  one-dimensional 
quadrature  rule.  The  error  in  integrating  /n(£)  over  (0,1)  with  these  points  is  less  than 
6  x  10-Cl.  The  method  has  been  shown  to  be  capable  to  handle  singularities  of  composite 
or  undetermined  types  [  Burton  1976  .  When  the  singularity  is  at  the  center  node  of  the  3 
noded  quadratic  element,  the  element  is  subdivided  such  that  the  quadrature  rule  can  be 
applied  on  either  side  of  the  node. 

In  three  dimensions  a  technique  given  by  ILachat  <£:  Watson  1976  was  used  to  integrate 
the  1/  r  singularity  in  Ut].  The  2x2  basis  square  in  -space  on  which  the  non-singular 
integrals  are  evaluated  is  subdivided  into  triangles,  the  singular  nodes  always  at  the  vertex. 
The  triangles  are  given  a  new  intrinsic  coordinate  system  (t?i.t72)  obtained  by  viewing  the 
triangle  as  a  degenerated  rectangle  in  the  (£i.sz)  space.  The  relationship  between  the  two 
sets  of  intrinsic  coordinates  is  given  in  terms  of  linear  shape  functions  defined  by 


Uv)  =  V  A',u,(^,1"1  M-7) 

a  —  1 


where  Na(t])  represent  the  linear  shape  functions.  These  triangular  subelements  in  the 
(r?  x ,  r/2 )  space  form  a  Jacobian  that  has  0(r)  behaviour.  The  0(  1  r)  singularity  of  the  kernel 
is  removed  numerically  when  multiplied  by  this  Jacobian  with  0(r)  behaviour. 


A. 3  Geometrical  Symmetry 

Symmetry  of  the  DAA  boundary  with  respect  to  coordinate  axes  is  accounted  for  within  the 
software.  This  is  implemented  by  reflecting  each  element  about  the  symmetry  axis  during 
the  construction  of  the  A  and  B  matrices  (equation  3.14).  Care  however,  is  required. when 
using  the  rigid  body  methodology  to  calculate  the  diagonal  terms  of  the  traction  kernel  TXJ. 
in  that  the  summation  of  the  off-diagonal  terms  must  be  performed  before  the  symmetry 
transformation  is  applied  to  each  component  of  TX].  Also,  the  displacements  and  tractions 
at  the  nodes  on  the  plane  of  symmetry  in  the  direction  across  the  plane  must  be  eliminated 
because  they  are  zero.  This  is  done  by  zeroing  the  corresponding  rows  and  columns  and 
by  placing  the  value  1.0  on  the  leading  diagonal. 
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SYMMETRIC  AND  NONSYMMETRIC 
MEDIUM  STIFFNESS  MATRICES 


The  accuracy  of  symmetric  and  nonsymmetric  medium  stiffness  matrices  is  evaluated  here 
by  computing  the  nodal  forces  generated  by  a  uniform  displacement  field  applied  to  a  spher¬ 
ical  cavity  in  an  infinite  elastic  medium.  The  correct  nodal  forces  follow  from  the  known 
traction  solution  Timoshenko  and  Goodier  1951  and  (3.16).  the  nodal  forces  produced  by 
the  nonsymmetric  stiffness  matrix  follow  from  (3.17).  and  the  nodal  forces  produced  by  the 
symmetric  stiffness  matrix  follow  from  (3.17)  with  Km  replaced  by  Km.  Figure  9  shows, 
for  the  discretization  of  Figure  5.  computational  error  in  nodal-force  magnitudes  computed 
with  the  symmetric  and  nonsymmetric  matrices:  Km  clearly  outperforms  Km.  It  should 
be  noted,  that  convergence  of  the  nodal  forces  generated  by  the  symmetric  medium  matrix 
Km  was  obtained  bv  successive  mesh  refinement. 
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FIELD  COMMAND  DEFENSE  NUCLEAR  AGENCY 
ATTN:  FCPR 
ATTN:  FCTKC  KELLER 
ATTN:  FCTT 
ATTN:  FCTTWSUMMA 
ATTN:  FCTXE 

JOINT  STRAT  TGT  PLANNING  STAFF 
ATTN:  JLK  (ATTN:  DNA  REP) 

ATTN:  JLKS 
ATTN:  JPST 
ATTN:  JPTM 

UNDER  SECYS  OF  DEF  FOR  RSCH  A  ENGRG 
ATTN:  STRATA  SPACE  SYS(OS) 

ATTN:  STRAT  A  THEATER  NUC  FOR  FVAJDA 

DEPARTMENT  OF  THE  ARMY 

BMD  ADVANCED  TECHNOLOGY  CENTER 
ATTN:  ATC-T 
ATTN:  ICRDAHB-X 

HARRY  DIAMOND  LABORATORIES 

ATTN:  DELHD-TA-L(81100)(TECH  LIB) 
ATTN:  SCHLD-NW-P 

U  S  ARMY  BALLISTIC  RESEARCH  LAB 
ATTN:  SLCBR  SS-T  (TECH  LIB) 


U  S  ARMY  CONCEPTS  ANALYSIS  AGENCY 
ATTN:  CSSA-ADL  (TECH  LIB) 

U  S  ARMY  CORPS  OF  ENGINEERS 
ATTN:  DAEN-ECE-T 
ATTN:  DAEN-RDL 

U  S  ARMY  ENGINEER  CTR  &  FT  BELVOIR 
ATTN:  DT-LRC 

U  S  ARMY  ENGINEER  DIV  HUNTSVILLE 
ATTN:  HNDED-SR 

U  S  ARMY  ENGINEER  DIV  OHIO  RIVER 
ATTN:  ORDAS-L  (TECH  LIB) 

U  S  ARMY  ENGR  WATERWAYS  EXPER  STATION 
ATTN:  J  STRANGE 
ATTN:  J  ZELASKO 
ATTN:  LIBRARY 

2  CYS  ATTN:  WESSD  J  JACKSON 
ATTN:  WESSE 

U  S  ARMY  MATERIAL  COMMAND 

ATTN:  DRXAM-TL  (TECH  LIB) 

U  S  ARMY  MATERIAL  TECHNOLOGY  LABORATORY 
ATTN:  TECHNICAL  LIBRARY 

U  S  ARMY  NUCLEAR  &  CHEMICAL  AGENCY 
ATTN:  LIBRARY 

USA  MISSILE  COMMAND 

ATTN:  DOCUMENTS  SECTION 

DEPARTMENT  OF  THE  NAVY 

LEAHY  (CG  16) 

ATTN:  WEAPONS  OFFICER 

NAVAL  FACILITIES  ENGINEERING  COMMAND 
ATTN:  CODE  04B 

NAVAL  POSTGRADUATE  SCHOOL 

ATTN:  CODE  1424  LIBRARY 

NAVAL  RESEARCH  LABORATORY 

ATTN:  CODE  2627  (TECH  LIB) 

NAVAL  SURFACE  WEAPONS  CENTER 
ATTN:  CODE  F31 

NAVAL  SURFACE  WEAPONS  CENTER 

ATTN:  TECH  LIBRARY  &  INFO  SVCS  BR 

NAVAL  UNDERWATER  SYSTEMS  CTR 

ATTN:  CODE  EM  J  KALINOWSKI 


DEPARTMENT  OF  THE  NAVY  (CONTINUED) 


STRATEGIC  AIR  COMMAND 
ATTN:  NRI/STINFO 


A 
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NAVAL  WAR  COLLEGE 

ATTN:  CODE  E  ll  (TECH  SERVICE) 

NAVAL  WEAPONS  EVALUATION  FACILITY 
ATTN:  CLASSIFIED  LIBRARY 

OFC  OF  THE  OEPUTY  CHIEF  OF  NAVAL  OPS 
ATTN:  N0P03EG 
ATTN:  NOP  981 

OFFICE  OF  NAVAL  RESEARCH 

ATTN:  CODE  474  N  PERRONE 

SPACE  &  NAVAL  WARFARE  SYSTEMS  CMD 
ATTN:  PME  117-21 

STRATEGIC  SYSTEMS  PROGRAMS  (PM-1) 
ATTN:  NSP-43  (TECH  LIB) 

DEPARTMENT  OF  THE  AIR  FORCE 

AFIS/INT 

ATTN:  I  NT 

AIR  FORCE  GEOPHYSICS  LABORATORY 
ATTN:  LWH/H  OSSING 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 
ATTN:  LIBRARY 

AIR  FORCE  OFFICE  OF  SCIENTIFIC  RSCH 
ATTN:  J  ALLEN 
ATTN:  W  BEST 

AIR  FORCE  SYSTEMS  COMMAND 
ATTN:  DLW 

AIR  FORCE  WEAPONS  LABORATORY.  AFSC 
ATTN:  NTE  M  PLAMONDON 
ATTN:  NTEDJ  THOMAS 
ATTN:  NTEDRHENNY 
ATTN:  SUL 

AIR  UNIVERSITY  LIBRARY 
ATTN:  AUL-LSE 

BALLISTIC  MISSILE  OFFICE/DAA 
ATTN:  EN 

ATTN:  MGEN  A  SCHENKER  ' 
ATTN:  MGENEFURBEE 
ATTN:  PP 

DEPUTY  CHIEF  OF  STAFF 
ATTN:  LEEE 

DEPUTY  CHIEF  OF  STAFF/AF-RDQI 
ATTN:  AF/RDQI 

FOREIGN  TECHNOLOGY  DIVISION.  AFSC 
ATTN:  NIIS  LIBRARY 

ROME  AIR  DEVELOPMENT  CENTER.  AFSC 
ATTN:  TSLD 

STRATEGIC  AIR  COMMAND 
ATTN:  INA 


STRATEGIC  AIR  COMMAND 
ATTN:  XPFS 

STRATEGIC  AIR  COMMAND 
ATTN:  XPQ 

DEPARTMENT  OF  ENERGY 

DEPARTMENT  OF  ENERGY 

ALBUQUERQUE  OPERATIONS  OFFICE 
ATTN:  CTID 
ATTN:  R  JONES 

DEPARTMENT  OF  ENERGY 

NEVADA  OPERATIONS  OFFICE 

ATTN:  DOC  CON  FOR  TECHNICAL  LIBRARY 

LAWRENCE  LIVERMORE  NATIONAL  LAB 
ATTN:  L-221D  GLENN 
ATTN:  L-53  TECH  INFO  DEPT  LIBRARY 

LOS  ALAMOS  NATIONAL  LABORATORY 

ATTN:  MS  P364  REPORTS  LIBRARY 

OAK  RIDGE  NATIONAL  LABORATORY 

ATTN:  CENTRAL  RSCH  LIBRARY 
ATTN:  CIVIL  DEF  RES  PROJ 

SANDIA  NATIONAL  LABORATORIES 

ATTN:  LIBRARY  &  SECURITY  CLASSIFICATION 
DIV 

SANDIA  NATIONAL  LABORATORIES 
ATTN:  ORG  7111  L  HILL 
ATTN:  ORG  7112  A  CHABAI 
ATTN:  TECH  LIB  3141 

OTHER  GOVERNMENT 

CENTRAL  INTELLIGENCE  AGENCY 
ATTN:  OSWR/NED 

DEPARTMENT  OF  THE  INTERIOR 
ATTN:  TECH  LIB 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

AEROSPACE  CORP 

ATTN:  LIBRARY  ACQUISITION  Ml/199 

APPLIED  RESEARCH  ASSOCIATES.  INC 
ATTN:  N  HIGGINS 

APPLIED  RESEARCH  ASSOCIATES.  INC 
ATTN:  SBLOUIN 

APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN:  DPIEPENBURG 

APPLIED  RESEARCH  ASSOCIATES.  INC 
ATTN:  R  FRANK 

AVCO  SYSTEMS  DIVISION 

ATTN:  LIBRARY  A830 

BDM  CORP 

ATTN:  AVITELLO 
ATTN:  CORPORATE  LIB 
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DEPT  OF  DEFENSE  CONTRACTORS  (CONTINUED) 

BDM  CORP 

ATTN:  F  LEECH 

BOEING  CO 

ATTN.  M/S  8K-22  D  CHOATE 

CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
ATTN:  T  AHRENS 

CALIFORNIA  RESEARCH  4  TECHNOLOGY,  INC 
ATTN:  K  KREYENHAGEN 
ATTN:  LIBRARY 
ATTN:  M  ROSENBLATT 
ATTN:  S  SCHUSTER 

CALIFORNIA  RESEARCH  4  TECHNOLOGY,  INC 
ATTN:  F  SAUER 

CALSPAN  CORP 

ATTN:  LIBRARY 

CARPENTER  RESEARCt  CORP 
ATTN:  HJ  CARPENTER 

UNIVERSITY  OF  DENVER 

ATTN:  J  WISOTSKI 

I  IT  RESEARCH  INSTITUTE 

ATTN:  DOCUMENTS  LIBRARY 
ATTN:  M  JOHNSON 

INSTITUTE  FOR  DEFENSE  ANALYSES 
ATTN:  CLASSIFIED  LIBRARY 

KAMAN  SCIENCES  CORP 
ATTN:  LMENTE 
ATTN:  LIBRARY 

KAMAN  SCIENCES  CORP 
ATTN:  LIBRARY 

KAMAN  TEMPO 

ATTN:  DASIAC 

KAMAN  TEMPO 

ATTN:  DASIAC 

LOCKHEED  MISSILES  4  SPACE  CO,  INC 

2  CYS  ATTN:  I  MATHEWS 
ATTN:  J  BONIN 

2 CYS  ATTN:  T GEERS 

LOCKHEED  MISSILES  4  SPACE  CO,  INC 
ATTN:  J  WEISNER 

MAXWELL  LABORATORIES,  INC 
ATTN:  J  MURPHY 

MCDONNELL  DOUGLAS  CORP 
ATTN:  RHALPRIN 

MERRITT  CASES,  INC 

ATTN:  J  MERRITT 
ATTN:  LIBRARY 

NEW  MEXICO  ENGINEERING  RESEARCH  INSTITUTE 
ATTN-  N  BAUM 


PACIFIC-SIERRA  RESEARCH  CORP 

ATTN:  H  BRODE.  CHAIRMAN  SAGE 

PATEL  ENTERPRISES,  INC 
ATTN:  M  PATEL 

R  4  D  ASSOCIATES 

ATTN:  C  K  B  LEE 
ATTN:  D  SIMONS 
ATTN:  J  LEWIS 
ATTN:  PHAAS 

ATTN:  TECHNICAL  INFORMATION  CENTER 
ATTN:  W  WRIGHT 

R  4  D  ASSOCIATES 

ATTN:  GGANONG 

RAND  CORP 

ATTN:  P  DAVIS 

RAND  CORP 

ATTN:  B  BENNETT 

S-CUBED 

ATTN:  DGRINE 
ATTN:  LIBRARY 
ATTN:  TRINEY- 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  TECHNICAL  LIBRARY 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  D  MAXWELL 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  WLAYSON 

SOUTHWEST  RESEARCH  INSTITUTE 
ATTN:  A  WENZEL 

SRI  INTERNATIONAL 

ATTN:  DKEOUGH 

STRUCTURAL  MECHANICS  ASSOC.  INC 
ATTN:  R  KENNEDY 

TELEDYNE  BROWN  ENGINEERING 
ATTN:  D  ORMOND 
ATTN:  F  LEOPARD 

TERRA  TEK,  INC 

ATTN:  S  GREEN 

TRW  ELECTRONICS  4  DEFENSE  SECTOR 
2  CYS  ATTN:  N  LIPNER 
ATTN:  PBHUTA 

ATTN:  TECHNICAL  INFORMATION  CENTER 

TRW  ELECTRONICS  4  DEFENSE  SECTOR 
ATTN:  E  WONG 
ATTN:  P  DAI 

WEIDLINGER  ASSOC,  CONSULTING  ENGRG 
ATTN:  TDEEVY 

WEIDLINGER  ASSOC,  CONSULTING  ENGRG 
ATTN:  M  BARON 

WEIDLINGER  ASSOC,  CONSULTING  ENGRG 
ATTN:  J  ISENBERG 
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